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Summary 

The performance of a three-dimensional Navier- 
Stokes solution technique in predicting the transonic 
flow past a nonaxisymmetric nozzle has been inves- 
tigated. The investigation was conducted at free- 
stream Mach numbers ranging from 0.60 to 0.94 and 
an angle of attack of 0°. Wind-tunnel data with the 
jet exhaust simulated by high pressure air were also 
obtained to compare with the numerical calculations. 
In the numerical procedure, the jet exhaust is repre- 
sented by a solid sting. 

‘ The numerical solution procedure employs 
the three-dimensional, time-dependent, Reynolds- 
averaged Navier-Stokes equations written in strong 
conservation form, a thin layer assumption, and the 
Baldwin-Lomax turbulence model. The equations 
are solved by using the finite-volume principle in 
conjunction with an approximately factored upwind- 
biased numerical algorithm. 

The Navier-Stokes solution technique predicts the 
flow over most of the afterbody quite well when no 
separation is present. Under these conditions, it does 
a particularly good job of capturing the negative 
pressure peak as the flow expands around the boattail 
shoulder. It also successfully predicts vortices which 
have been experimentally shown to trail from the 
nonaxisymmetric nozzle. As expected, the procedure 
does not predict the flow well at the juncture of the 
nozzle exit and the solid jet plume simulator. This 
result reiterates the well-established fact that a solid 
plume simulator does not adequately represent the 
jet exhaust. With the Baldwin-Lomax turbulence 
model, the numerical solution procedure also gives 
poor results when there are strong shocks and large 
separated areas on the nozzle. 

Introduction 

When a typical fighter airplane is flying at tran- 
sonic speeds, 40 to 50 percent of its drag is associated 
with its afterbody (refs. 1 and 2). In the afterbody 
region, flows over the fuselage and empennage merge 
and further interact with the propulsive exhaust to 
create a complex flow field which strongly influences 
the aerodynamic characteristics of the airplane. The 
flow field is rotational and highly interactive; hence 
viscous effects are very important, especially at off- 
design conditions. 

Traditionally, engineers have utilized wind tun- 
nels to investigate this complicated flow field and de- 
rive configurations with desirable flight characteris- 
tics. Theoretical methods which would adequately 
handle the intricate flow phenomena did not ex- 
ist. Now, however, numerical solution techniques 
which will compute complicated aerodynamic flows 


are becoming available and are being used to com- 
plement, but not replace, wind tunnels as inves- 
tigative tools. In general, computational methods 
more readily yield information about the complete 
flow field and isolated components than wind tun- 
nels; thus, they are well suited for use in gaining an 
understanding of the mechanics of the flow and iden- 
tifying components with undesirable aerodynamic 
characteristics. 

Inviscid, irrotational computational methods 
have already proven to be useful tools in the design 
of wings at cruise conditions. (See refs. 3 and 4.) 
However, for predicting the complex flow field gen- 
erated in the afterbody region at transonic speeds, 
Navier-Stokes techniques are needed to adequately 
simulate the physics of the flow. Examples of Navier- 
Stokes calculations which have been made for af- 
terbody flows include those of Swanson and Deiw- 
ert. Swanson solved the axisymmetric form of the 
equations for the flow over axisymmetric afterbodies 
(ref. 5), and the two-dimensional form of the equa- 
tions for the flow over nonaxisymmetric afterbodies 
(ref. 6). He solved the three-dimensional equations 
for the internal flow only (ref. 6). These compu- 
tations were made on a CDC® CYBER 203 com- 
puter using an alternating-direction-implicit numer- 
ical algorithm with central differencing to integrate 
the equations over time and space. 

Deiwert solved the axisymmetric form of the 
Navier-Stokes equations for both transonic and 
supersonic flow over axisymmetric afterbodies (refs. 7 
and 8). He solved the three-dimensional form of 
the equations for the supersonic flow past a body 
of revolution with a centered propulsive jet at angle 
of attack (ref. 9). These solutions were obtained on 
a CDC® CYBER 205 computer. Deiwert also used 
an alternating-direction-implicit numerical algorithm 
with central differencing to integrate the equations. 

Recently, advances such as upwind differencing 
have been made in numerical algorithms for solv- 
ing the three-dimensional Navier-Stokes equations. 
See references 10 through 13, for instance. In refer- 
ence 13, an upwind-biased numerical algorithm was 
used to solve these equations for the flow over a pro- 
late spheriod at low subsonic Mach numbers and 
at angle of attack. Although the emphasis was on 
the whole spheroid, afterbody results were presented. 
However, in none of these investigations were the 
three-dimensional Navier-Stokes equations solved for 
the external afterbody flow at transonic speeds. 

The present paper examines the performance of 
this upwind-biased numerical algorithm in solving 
the three-dimensional Navier-Stokes equations for 
the subsonic and transonic flow over a nonaxisym- 
metric nozzle. It extends the previous work in three 


main ways. First, the three-dimensional Navier- 
Stokes equations are solved for the flow past an 
afterbody at transonic speeds. Second, the upwind- 
ing algorithm is applied to this specific problem. Fi- 
nally, the application is for a nonaxisymmetric noz- 
zle typical of those being advocated for advanced 
fighter airplanes. The free-stream test Mach num- 
bers ranged from 0.60 to 0.94, and the angle of at- 
tack was 0°. Wind-tunnel data to compare with the 
numerical calculations were also obtained. 

The Navier-Stokes solution procedure used in 
making the calculations is described in references 10 
through 13. In the procedure, the three-dimensional, 
time-dependent, Reynolds- averaged Navier-Stokes 
equations are written in strong conservation form and 
a thin layer assumption is made to simplify the dis- 
sipation terms. For the case of a turbulent bound- 
ary layer, these terms are simulated w T ith a Baldwin- 
Lomax turbulence model. The solution procedure 
operates on the finite-volume principle and employs 
an approximately factored upwind-biased numerical 
algorithm to integrate the equations over time and 
space. For the present calculations, the jet exhaust 
is represented by a solid sting. 

The experimental portion of the investigation was 
conducted in the Langley 16-Foot Transonic Wind 
Tunnel. The data presented are part of a broad 
experimental data base which was obtained for the 
nonaxisymmetric variable-geometry nozzle described 
previously. In the wind-tunnel experiment, high 
pressure air was used to simulate the jet exhaust. 

The performance of the numerical Navier-Stokes 
procedure in predicting the three-dimensional after- 
body flow field is assessed over the range of condi- 
tions investigated. An assessment is made of the 
suitability of the Baldwin-Lomax turbulence model 
to handle the turbulent diffusion in the flow regimes 
encountered. 
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Cf 

Cl 

Cd 


Jacobian matrix of transformed flux 
vector, F 

reference area, A re f = 1.0 
dimensions (see fig. 4(c)) 
skin friction coefficient, 

Yoo 

lift coefficient, -j 
drag coefficient, 


C p pressure coefficient, p 

Cp critical pressure coefficient 
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F, G, H 


F 

G 

H 

H 

H v 

I 

i,j,k 

i,j,k 
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M 0 
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P 

Q 
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Q 

q 

Qoc 


R 


speed of sound 

specific heat at constant pressure 
drag 

total energy per unit volume 

flux vector in physical (Cartesian) 
coordinate direction x, y, and z, 
respectively 

flux vector in £ transformed coordi- 
nate direction 

flux vector in rj transformed coordi- 
nate direction 

total enthalpy per unit mass 

flux vector in £ transformed coordi- 
nate direction 

viscous flux vector in £ transformed 
coordinate direction 

identity matrix 

index in £, y, and £ transformed 
coordinate direction, respectively 

unit vector in x,y, and z physical 
coordinate direction, respectively 

Jacobian of transformation from 
physical coordinate system to 
computational coordinate system 

lift 

reference length, length of model 
from nose to end of solid jet plume 
simulator 

length of model from nose to 
juncture of nozzle exit and solid jet 
plume simulator 

free-stream Mach number 

number of time steps 

pressure 

vector of dependent flow variables 

transformed vector of dependent 
flow variables 

square of velocity (see eqs. (3)) 
vector of primitive variables 
free-stream dynamic pressure, 

^Po o{ u ce v oc 4 " w oc) 

residual vector 
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Rl Reynolds number based on length of 

model from nose to nozzle exit 

S directed surface area of cell face 

T diagonalizing matrix 

t time 

U contravariant velocity component in 

£ transformed direction 

u normalized x velocity component in 

physical coordinate system 

V contravariant velocity component in 

r 7 transformed direction 

v normalized y velocity component in 

physical coordinate system 

v volume 

W contravariant velocity component in 

£ transformed direction 

w normalized z velocity component in 

physical coordinate system 

x axial physical (Cartesian) coordi- 

nate, origin at nose of model 

y horizontal physical (Cartesian) 

coordinate, origin at nose of model 

z vertical physical (Cartesian) coordi- 

nate, origin at nose of model 

P exponent (see fig. 4(c)) 

7 ratio of specific heats 

A incremental quantity or forward 

difference formula 

A t incremental time 

, Sfj , differencing operator in £, 77, and 

£ direction, respectively (see for 
example eq. (13)) 

£ radial coordinate in transformed 

coordinate system 

77 circumferential coordinate in 

transformed coordinate system 

9 circumferential coordinate, radians 

k spatial differencing parameter 

(eq. (23)) 

K e effective conductivity 

li viscosity 


€ 

axial coordinate in transformed 
coordinate system 

p 

density 

a 

laminar Prandtl number 

<?t 

turbulent Prandtl number 


local skin friction at wall 

V 

Subscripts: 

backward difference formula 

c 

corner dimension (fig. 4(c)) 

e 

effective 

i, j, k 

index of cell location in £, 77, and £ 
transformed coordinate direction, 
respectively 

L 

left 

n 

nozzle (see fig. 4) 

R 

right 

t 

turbulent 

00 

free-stream conditions 


A boldface symbol denotes a vector quantity or 
a matrix; a tilde (~) over a symbol denotes a Roe- 
averaged quantity; a caret C) over a symbol denotes 
a transformed quantity; the superscript T denotes a 
transposed matrix. 

Unless otherwise noted, all variables are non- 
dimensionalized by appropriate combinations of the 
free-stream parameters and a reference length of 
1 inch. 

Computational Procedure 

The three-dimensional Navier-Stokes numerical 
solution procedure used is described in references 10 
through 13. In the procedure, the three-dimensional, 
time-dependent, Reynolds- averaged Navier-Stokes 
equations are written in strong conservation form. 
A thin layer assumption is made to simplify the dis- 
sipation terms. The turbulent contribution to these 
terms is simulated with the Baldwin-Lomax turbu- 
lence model (ref. 14). 

The solution procedure employs the finite-volume 
principle where the spatial derivatives in the equa- 
tions are evaluated as conservative flux bfJances 
across the grid cells. The fluxes at the cell interfaces 
are determined with upwind-biased flux-difference 
splitting in combination with a gradient limiting 
procedure to ensure monotonicity across disconti- 
nuities such as shock waves. The scheme is spa- 
tially third-order accurate. The time differencing al- 
gorithm used in the computational procedure is an 
approximately factored alternating-direction-implicit 
scheme in delta form. For completeness, a detailed 
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description of the computational procedure, includ- 
ing equations, used in this paper is presented in the 
appendix. 

Computational Domain and Grid 

A sketch of the three-dimensional computational 
domain illustrating the domain boundaries, grid 
topology, and model is presented in figure 1. To min- 
imize computer run time and storage requirements, 
symmetry is assumed about the vertical plane, and 
only one half of the cylindrical domain is computed. 
The computational grid is body fitted (grid lines co- 
incide with nacelle surface and other boundaries) in 
order to facilitate implementation of the boundary 
conditions. 

The transverse, or cross-flow, planes of the three- 
dimensional coordinate system are generated by us- 
ing the Thompson-Thames-Mastin elliptic method of 
reference 15. To construct the three-dimensional sys- 
tem, a series of these two-dimensional grids are posi- 
tioned normal to an axis which is nominally aligned 
with the free-stream flow direction. This method of 
construction results in a sheared, H-type grid topol- 
ogy in the longitudinal planes, and an O-type topol- 
ogy in the cross-flow planes. (See fig. 1.) 

The spacing of the mesh is stretched away from 
the model, with grid lines being radially clustered 
near the body surface and axially clustered on the 
afterbody. The spacing normal to the model sur- 
face is illustrated in figure 1. Figure 2 presents a 
wire- frame sketch of the model showing a typical grid 
spacing on its surface for a coarse mesh. The fore- 
body, afterbody, and plume simulator are also de- 
fined in the figure. 

Calculations were made with two computational 
meshes of different densities. The first grid system 
consisted of 64 grid planes in the axial direction (15 
on the afterbody), 32 in the radial direction, and 10 
in the circumferential direction. To assess the effects 
of grid refinement, calculations were also made with 
a mesh having 129 planes in the axial direction (29 
on the afterbody), 65 in the radial direction, and 33 
in the circumferential direction. For both meshes, 
the first radial grid point off the model surface is 
located 0.00076 inch from the surface. Based on 
the length of the model from the nose to the nozzle 
exit (approximately 63.25 inches), this amounts to 
0.000012 of the model length. Most of the results 
presented were computed on the finer mesh. 

Experimental Apparatus and Procedure 

The experimental portion of this investigation 
was conducted in the Langley 16-Foot Transonic 
Wind Tunnel which is a continuous, atmospheric, 


single-return wind tunnel with an octagonal, slotted 
test section. A thorough descripton of the wind 
tunnel can be found in reference 16. The data 
presented are part of a broad experimental data base 
which was obtained for a nonaxisymmetric variable- 
geometry nozzle. The nozzle is typical of those 
currently advocated for advanced fighter airplanes. 
The parts of the data base that have already been 
published can be found in references 17 and 18. 

For testing, the nozzle was mounted on a generic 
forebody which was supported in the wind tunnel 
by a sting-strut arrangement. Figure 3 shows a 
photograph of the model mounted in the wind tunnel. 
Details of the forebody and mounting system can 
be found in reference 17. Surface pressures on the 
nozzle were measured at free-stream Mach numbers 
up to 1.3. High pressure air was used to simulate 
the jet exhaust over a range of nozzle pressure ratios. 
During the investigation, Reynolds numbers based on 
the overall model length of 63.3 inches ranged from 

17.2 x 10 6 at a free-stream Mach number of 0.60 to 

29.3 x 10 6 at a free-stream Mach number of 0.94. 

The specific nozzle configuration considered in 

this investigation represents a transonic-cruise, dry- 
power setting of the variable- geometry nozzle. Fig- 
ure 4 presents a sketch of the configuration and 
gives its pertinent dimensions. Externally, the con- 
figuration is 7.24 inches long, 6.20 inches high, and 
6.80 inches wide. At each longitudinal station, it 
has straight external sides with superelliptic corners. 
Both the external sidewall centerline and flap cen- 
terline are composed of circular arcs tangent to a 
straight section (refer to figs. 2, 3, and 4), with the 
circular-arc radii equal to 10.86 inches. The side- 
wall boattail angle is 6.93°, and the flap boattail an- 
gle is 17.56°. Internally, the sidewalls are flat. The 
nozzle has a throat area of 2.57 inch 2 , a throat as- 
pect ratio of 2.38, and an exit aspect ratio of 1.91. 
Its expansion ratio is 1.25 which yields a total pres- 
sure ratio for perfect expansion of approximately 
4, the pressure ratio used for comparison with the 
computations. 

Results and Discussion 

Navier-Stokes solutions were obtained for the 
nonaxisymmetric afterbody at free-stream Mach 
numbers of 0.60, 0.80, and 0.94 and an angle of at- 
tack of 0°. To obtain the solutions, the afterbody 
was coupled with a forebody; thus, a complete model 
was created which matched the wind-tunnel configu- 
ration. The Reynolds numbers corresponding to the 
three free-stream Mach numbers were, respectively, 
17.2 x 10 6 , 19.5 x 10 6 , and 29.3 x 10 6 based on the 
model length. All calculations were made with the 
three-dimensional Navier-Stokes procedure described 
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previously. In the experimental portion of the in- 
vestigation, pressures were measured on the model 
afterbody. The computed pressure coefficients are 
compared with the wind-tunnel data. 

Validity of the Solutions 

Numerical convergence . In order to help assess 
the validity of the solutions, the numerical conver- 
gence characteristics were monitored as the solutions 
progressed in time. Figures 5 and 6 present the con- 
vergence history and the general features of the flow 
for solutions at free-stream Mach numbers of 0.60 
and 0.80. At a Mach number of 0.60 (fig. 5(a)), the 
residual drops approximately four orders of magni- 
tude in 1000 time steps. The drag coefficient reaches 
a constant value in approximately 400 time steps. 
The lift coefficient is zero throughout its time his- 
tory since the computations are made at an angle 
of attack of 0°. These three histories indicate that 
the solution is converging and that it has essentially 
reached a steady state in about 400 time steps. 

Figure 5(b) shows the computed pressure coeffi- 
cient distributions along the top and side of the com- 
plete body. Comparing them with wind-tunnel data 
included in the figure shows that the calculated pres- 
sures along the entire length of the nacelle surface are 
reasonable. They agree well with the experimental 
data except at the compression peak near the junc- 
ture of the afterbody and the solid jet plume simula- 
tor. This discrepancy is discussed more fully later. 

Figure 6 presents the convergence history and the 
solution at a Mach number of 0.80. The residual 
drops several orders of magnitude to a plateau where 
it oscillates with a small amplitude about a constant 
value. The plateau effect is thought to be a result 
of the discontinuous derivative of the min mod flux 
limiter in certain regions of the flow (ref. 11). There 
is good correlation between the computed pressures 
and wind-tunnel data over most of the nacelle. Again 
the exception to the good correlation is at the junc- 
ture of the afterbody and the plume simulator. 

Grid convergence . A brief study of the effect 
of grid resolution on the quality of the solution was 
made in order to further assess the validity of the so- 
lution. In making this study, calculations were made 
with two computational meshes of different densities. 
The coarser mesh consisted of 64 grid planes in the 
axial direction (15 on the afterbody), 32 in the ra- 
dial direction, and 10 in the circumferential direction. 
The finer mesh had 129 planes in the axial direction 
(29 on the afterbody), 65 in the radial direction, and 
33 in the circumferential direction. Figure 7 presents 
the pressure distributions along the top and side of 


the afterbody and the cross-plane pressure distribu- 
tions at the nozzle exit for both calculations. Both 
solutions agree well with the experimental data on 
the top and side rows except near the exit (figs. 7(a) 
and 7(b)). Although using the finer grid results in 
slightly better resolution for the top row, a compar- 
ison of the two solutions shows that grid resolution 
is not an issue on the top and side rows. The cross- 
plane pressure distributions at the exit (presented in 
fig. 7(c)), however, show that a fine grid is required 
if the solution is to be adequately resolved at the 
corners of the nonaxisymmetric nozzle. 

The residual histories presented at free-stream 
Mach numbers of 0.60 and 0.80 (figs. 5 and 6) in- 
dicate that the solutions are converging numerically. 
The pressures at these Mach numbers look reason- 
able in all regions of the flow and agree well with 
the wind-tunnel values in regions where experimen- 
tal data were obtained except near the exit. Fur- 
thermore, the solutions appear to be grid converged 
for the vertical and horizontal planes containing the 
model axis. Thus, within the constraints of the math- 
ematical modeling of the physics, the solutions in the 
vertical and horizontal longitudinal planes seem to be 
converged to reasonably accurate representations of 
the true time-averaged flows being calculated. 

Overall Solution Features and Comparison 
With Experiment 

Figures 8 through 10 present details of the solu- 
tions for the three free-stream test Mach numbers. 
The details are specifically aimed at the afterbody 
region. The figures present Mach number contours, 
pressure distributions, and skin friction distributions 
for the top and side longitudinal planes or rows as the 
case may be. The Mach number contours are shown 
in parts (a) and (b) of the figures, the pressures in 
parts (c) and (d), and the skin friction distributions 
in parts (e) and (f). 

Moo = 0.60. At a free-stream Mach number of 
0.60, presented in figure 8, the flow is relatively mild 
with no shocks or separated regions present. The 
general structure of the flow field can be seen by 
looking at the local Mach number contours for the 
top and side longitudinal planes, shown in parts (a) 
and (b) of the figure. In particular, one can see the 
thinning of the boundary layer where the flow ac- 
celerates around the shoulder of the boattail and its 
corresponding thickening at the nozzle exit where the 
boattail ends and the solid jet plume simulator be- 
gins. The boundary- layer effects are more prominent 
along the top of the nozzle than they are on the side 
of the nozzle. This characteristic of the flow is to be 
expected since the boattail angle is much greater on 
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the top and bottom of the nozzle than on its sides. 
The larger angles result in greater accelerations of 
the flow accompanied by more severe changes in the 
boundary layer. 

The pressure distributions shown in figures 8(c) 
and (d) reflect the structure of the flow indicated by 
the Mach number contours. They also correlate well 
with the wind-tunnel data over most of the nozzle 
boattail. A comparison of the calculated pressures 
and the data shows that the computations capture 
the level of the suction pressure peak exceptionally 
well. This good correlation occurs on both the top 
and side of the nozzle. 

The one exception to the good correlation is at 
the juncture of the nozzle exit and the solid jet 
exhaust plume simulator. The discrepancy at the 
juncture is thought to be predominantly caused by 
the solid plume simulator inadequately modeling the 
jet exhaust plume. These results reiterate the well- 
established fact that the jet exhaust must be cor- 
rectly modeled in order to accurately predict the flow 
near the nozzle exit. (See refs. 19 through 23, for ex- 
ample.) Not only does the solid plume simulator in- 
accurately model the shape of the jet exhaust plume, 
but it also fails to account for any entrainment of the 
afterbody flow by the jet exhaust. 

The importance of accurately simulating both the 
blockage effects of the jet exhaust plume on exter- 
nal flow and the entrainment produced by the jet 
exhaust was experimentally demonstrated in refer- 
ences 19 and 20. In reference 21, entrainment was 
accounted for with some success when using solid jet 
plume simulators by reducing the size of the sim- 
ulator. However, the size of the simulator needed 
for best results depended on the flow conditions. 
In references 22 and 23, the viscous-inviscid inter- 
action computational procedure developed to simu- 
late jet entrainment gave good predictions for flows 
with weak interactions. However, in the general case, 
the most accurate way to account for all the effects 
of the jet exhaust would be to include the exhaust 
in the Navier-Stokes calculations as was done in 
reference 9. 

Another factor which may contribute to the calcu- 
lated pressures being higher than the data is the in- 
ability of the turbulence model to predict the dissipa- 
tion in the region of severe curvature at the juncture 
of the nozzle and simulator. (See refs. 24 and 25.) 
This factor is considered minor in this case though. 

The calculated skin friction distributions are pre- 
sented in figures 8(e) and (f). No experimental skin 
friction data exist for the afterbody. The skin fric- 
tion levels seem to be reasonable and correlate with 
the pressure distributions. 


Moo = 0.80. The solution at a free-stream Mach 
number of 0.80 is presented in figure 9. At this Mach 
number, the flow field has the same general structure 
as it does at a Mach number of 0.60, but the gradients 
and general features of the flow are magnified. A look 
at the pressure distributions presented in figure 9(c) 
shows that at a free-stream Mach number of 0.80, the 
flow goes slightly supersonic as it expands around the 
shoulder on top of the afterbody. Even under these 
conditions, the computations capture the negative 
pressure peak at the boattail shoulder exceptionally 
well for both the top and side of the nozzle. The 
skin friction predictions look reasonable at this Mach 
number also. 

Again, as expected, the flow is not predicted 
well at the juncture of the nozzle and the solid 
plume simulator. The same comments concerning 
the inaccurate simulation of the jet plume apply 
at this Mach number as at Mach 0.60. Hence, 
the results give strong motivation to include the jet 
exhaust and its interaction with the external flow in 
the Navier-Stokes calculations. 

Moo = 0.94 • The results of the computations 
at a free-stream Mach number of 0.94 are presented 
in figure 10. At this Mach number, the solution 
procedure predicts large supersonic regions of flow. 
A shock exists on the top of the nozzle boattail and 
either a very weak shock or a strong compression 
exists on its side as can be seen from parts (a) and 
(b) of the figure. A much greater thinning of the 
boundary layer can be seen as the flow accelerates 
around the boattail shoulder than can be seen at 
the two lower Mach numbers. In addition, the 
predicted thickening of the boundary layer as the flow 
recompresses is much more dramatic than was found 
at the two lower Mach numbers. 

Unfortunately, a comparison of the calculated 
pressures and the experimental data shows that the 
predicted flow over most of the top of the nozzle is 
only qualitatively correct in a very general sense. The 
measured wind-tunnel pressures show that the shock 
on top of the nozzle occurs much farther upstream 
than is predicted. Also, the data indicate that a large 
area on top of the nozzle aft of the shock is separated 
with a pressure plateau in the separated region. The 
calculated pressures contain no significant pressure 
plateau and, thus, by themselves indicate no separa- 
tion aft of the shock. The computed pressure recov- 
ery in this region is also much too large. 

Improper simulation of entrainment undoubtedly 
has some role in overpredicting the pressure recovery 
in this case. However, the predominant reason for 
the poor prediction of the flow at this Mach number 
is attributed to the use of an algebraic eddy viscosity 
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turbulence model. (See refs. 24 and 25.) Reynolds- 
averaged Navier-Stokes calculations using eddy vis- 
cosity turbulence models typically predict the shock- 
wave too far downstream in transonic flows. They 
also commonly fail to give a significant plateau in the 
separated region. Some researchers have used relax- 
ation equations for the eddy viscosity and achieved 
plateaus (for example, ref. 5). However, the results 
depended on the relaxation length, and the optimum 
length often depends on the flow. (See ref. 24 for 
a discussion of these points.) Hence, other classes 
of turbulence models need to be investigated for cal- 
culating transonic propulsion integration effects for 
fighter-type configurations. 

Flow-Field Vortices 

Formation and structure , Computed afterbody 
pressures along the top and side rows of the afterbody 
are compared in figure 11 for the turbulent solution at 
a free-stream Mach number of 0.80. In the expansion 
region, the pressures on top of the afterbody are 
considerably lower than the pressures on its side. 
Conversely, over the last part of the boattail, the 
recompression region, the pressures on the top are 
higher. 

The impact of the circumferential variation of the 
pressures on the overall flow pattern is indicated in 
figure 12. The figure presents computed particle 
traces for the grid cells immediately next to the 
afterbody surface. Particle traces such as these give 
a picture of the flow patterns just off the surface and 
hence can be considered the numerical equivalent of 
experimental oil flows. An examination of the traces 
shows that the initially higher pressures on the side 
of the afterbody cause the flow to sweep from the 
sides, around the corners, and toward the center of 
the top of the afterbody. As the flow approaches the 
exit, the reversal of the circumferential location of 
the high and low pressures changes the direction of 
the sweep of the flow. 

This swirling behavior of the flow suggests that 
vortices are likely to be generated. In order to 
determine if vortices exist, consider the cross-flow 
velocity vectors at several axial stations along the 
afterbody and downstream of the exit. Plots of these 
vectors are presented in figure 13. They show that at 
the beginning of the boattail, part (a) of the figure, 
the streamlines are mainly in the axial direction. 
Farther down the afterbody, strong inward radial 
velocities are created. Also, the flow begins to sweep 
from the sides of the afterbody, around the corners, 
and toward the centerline of the top and bottom of 
the afterbody in the manner indicated by the particle 


traces. At the nozzle exit, part (d), the suggestion of 
vortices at the corners of the nozzle can be seen. 

At stations aft of the exit, the vortices become 
well formed and appear to move slightly outward 
as the flow progresses downstream. Finally, about 
lj nozzle lengths aft of the exit, they appear to be 
getting larger in size and much less distinct which 
indicates that they are dissipating. In addition, at 
this point, the radial velocities have diminished in 
magnitude and changed direction from inward to 
outward. 

Thus, the numerical solution technique predicts 
the formation and demise of vortices in the vicin- 
ity of the nonaxisymmetric nozzle. The existence 
of vortices behind this nonaxisymmetric nozzle was 
also discovered in the experimental investigation of 
reference 18. 

Significance of vortices. The generation of vor- 
tices behind this nonaxisymmetric nozzle suggests 
that excessive drag is being created. Possibly, an 
afterbody with equal pressures on the sides and top 
and bottom where the flow does not tend to sweep 
around the corners would generate either weaker vor- 
tices or no vortices at all and hence have less drag. 
This idea is experimentally corroborated by the re- 
sults of reference 26. During that experimental inves- 
tigation, it was found that a nonaxisymmetric nozzle 
with equal boattail angles on the sides and top and 
bottom had less drag than the other nonaxisymmet- 
ric nozzles tested. 

Computer Resources Required 

The three-dimensional Navier-Stokes numerical 
solution technique takes approximately 75 micro- 
seconds per grid point per time step on the, 
CYBER VPS 32 computer. This amounted to ap- 
proximately 7.8 hours of computer time to get a 
steady-state solution on the refined grid. While this 
may seem like a large amount of computer time, the 
solutions yield a large amount of information about 
the complete flow field as has been demonstrated. 

Concluding Remarks 

The performance of a three-dimensional Navier- 
Stokes solution technique in predicting the transonic 
flow past a nonaxisymmetric nozzle has been inves- 
tigated. The investigation was conducted at free- 
stream Mach numbers ranging from 0.60 to 0.94 and 
an angle of attack of 0°. Wind-tunnel data with the 
jet exhaust simulated with high pressure air were also 
obtained to compare with the numerical calculations. 
In the numerical procedure, the jet exhaust is repre- 
sented by a solid sting. 
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The numerical solution procedure employs the 
three-dimensional, unsteady, Reynolds-averaged 
Navier-Stokes equations written in strong conserva- 
tion form, a thin layer assumption, and the Baldwin- 
Lomax turbulence model. The equations are solved 
by using the finite-volume principle in conjunction 
with an approximately factored upwind-biased nu- 
merical algorithm. 

The numerical experiments demonstrate that at 
the free-stream Mach numbers of 0.60 and 0.80, 
the Navier-Stokes solution technique predicts the 
flow over most of the afterbody quite well when 
no separation is present. At these Mach numbers, 
the solution technique does a particularly good job 
of capturing the negative pressure peak as the flow 
expands around the boattail shoulder. 

The Navier-Stokes technique also successfully pre- 
dicts the generation of vortices by the nonaxisymmet- 
ric afterbody, a phenomenon which has been exper- 
imentally verified. This successful prediction sheds 
insight into possible ways to reduce the afterbody 
drag and, thus, illustrates the potential useful- 
ness of numerical techniques in solving propulsion- 
integration problems. 


As expected, the procedure does not predict the 
flow well at the juncture of the nozzle exit and the 
solid jet plume simulator. This result reiterates 
the well-established fact that a solid plume simula- 
tor does not adequately represent the jet exhaust. 
Hence, it gives strong motivation to include the ex- 
haust in the Navier-Stokes calculations. The numer- 
ical procedure also does not predict the flow well at 
the free-stream Mach number of 0.94 where there are 
strong shocks and large separated areas on the noz- 
zle. Since algebraic eddy- viscosity turbulence mod- 
els generally fail for this category of fluid flows, the 
suitability of other classes of turbulence models for 
calculating transonic afterbody flows needs to be 
investigated. 

Nevertheless, considering all the results above, 
the three-dimensional Navier-Stokes numerical tech- 
nique in its present form is a very useful research tool. 
Combining the use of the solution technique with the 
use of the wind tunnel to investigate transonic after- 
body flows should yield new insights into the physics 
involved in these flows. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
April 28, 1989 
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Appendix 


Computational Procedure 


Governing Equations 

The basic mathematical model of the physical laws governing the motion of a compressible 
fluid with dissipation is the set of Navier-Stokes equations. In the solution procedure used 
in this investigation, the three-dimensional, time-dependent, Reynolds- averaged Navier-Stokes 
equations are normalized and written in strong conservation form for a Cartesian coordinate 
system (x, y,z). Body forces play an insignificant role in the afterbody flow problem and are 
neglected. Since the dominant dissipative effects for the problem arise mainly from diffusion 
normal to the afterbody surface, a thin-layer assumption is made by retaining only those 
diffusion terms which are normal to the body surface. The resulting time-dependent equations 
for conservation of mass, linear momentum, and energy can then be expressed in terms of a 
fixed generalized coordinate system (£,77, () as 


where 




r p 1 

pu 

pv 

pw 


F 


r pu \ 


1 


pUu + £ x p 
< pUV+tyP > 
pUw + £ z p 
. ( e+p)U j 



( py 

pVu + T) x p 
< pVV + T]yP > 
pVw + T] z p 

, ( e + p)V . 


(2a) 


(2b) 


(2c) 
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r pw \ 



pWu + C ; xP 
pWv + ( y p > 
pWw + C z p 
(e + p)W , 


(2d) 


0 


H„ = 


M<x> Pe 


l 01 (£) 


<Ai«c + Cx02 
01 «C + Cj/02 
01W< + Cz02 

- & ^Ti (c 2 ) + M, 'fc J 


£ c p Pe (t 


(2e) 


In these equations u , u, and w are velocity components in the physical coordinate system 
(x, ?/, > 2 ;), p is the density, p is the pressure, c p is the specific heat at constant pressure, and e is 
the total energy per unit volume. The terms p e and K e are, respectively, the effective viscosity 
and the effective conductivity which are composed of laminar and turbulent elements. Also, 


01 — C a; + Cy + C z 

02 = ^ (Cr^c + Cy v c + Cz w c) 

2 2,2, 2 

^ = u + V + W 




PooQoo^ 

Poo 


(3) 


The term £ is the reference length for determining the Reynolds number. In this paper, it is 
the length of the model from the nose to the exit (i.e., from the nose to the juncture of the 
nozzle and the solid jet plume simulator (fig. 2)). 

The contravariant velocity components used in equations (2) are defined as 
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U = ixU + tyV + £, Z W 
V = T] x u + T}yV + T) z w 
W = CxU + C,yV + C zW 


( 4 ) 



and the metrics of the transformation as 


Zx = J (VriH - VC, z v) ' 

vx = J(vch ~ Vt z > <) 

Zx = J(y^Zr, - y v z {) 

Zy = J ( x <; z y - x v z 0 

Vy = J(x^z (-x^Zf) > 

Z y = J(xt)Z£ ~ x^Zff) 

Zz = J(xr,y^ - x c y v ) 

Vz = J {x^ - x 
Z z = j {x^y n - x n yz) 

a 

where J is the Jacobian of the transformation and is given by 

J~ l = x^(y n z^ - y^Zjj) - y^x^ZQ - X(Z V ) + z^x^ - y^x^) 

Finally, the relations between the energy, pressure, and enthalpy for an ideal gas 


( 5 ) 


( 6 ) 


P = { 7- 1) (e- ^pg 2 ) | 

> (7) 

J 

complete the system of equations. 

In the momentum equations, the effective viscosity fi e is computed as 

Ve = M + fit — f 1 (8) 

where fi is the molecular viscosity and fit is the turbulent viscosity. Similarly, in the energy 
equation, the effective conductivity K e is given by 


Cp C p Cpfi 

«e = « + K t = —fi + —fit = — — 

a at a 



(9) 


where a is the laminar Prandtl number and at is the turbulent Prandtl number. The turbulence 
quantities are evaluated by using the algebraic eddy viscosity turbulence model of Baldwin and 
Lomax (ref. 14). 

The governing equations are written here in a finite-difference form to facilitate describing 
the numerical algorithm. However, the equations are solved with the finite- volume technique or 
by summing the fluxes crossing the faces of each volume formed by the grid. Any appropriate 
pressure-area terms acting on the faces are added. The result is then equated to the time 
rate of change of the conserved quantity in each respective volume. In the finite-volume sense, 


the state Qij,k ls n °t regarded as an approximation to the state at the cell center but as an 
approximation to the average state in the cell (ref. 27) or as 


i» = j_ r 

li ' j ' k Av J (i i 


(i+2)A^,(j+j)Ar;,(fc+2)AC 


(i-j)A^,(i-^)Ar?,(fc-^)AC 


Q(£,/?,C,nA<) dv 


(10) 


where Av is the incremental volume, At the incremental time, and n denotes the time step. In 
the present formulations, the spacing in the transformed plane is unity; that is, 


Z i+ = 1 l 


vrvi = 1 


C fc+ 1 -C fc _i = 1 


(11) 


The relationship between the finite-volume and finite- difference representations of the 
equations can be seen by comparing the discretized forms of each representation (ref. 28). 
Specifically, this comparison illustrates how the Jacobian and metric derivatives in the finite- 
difference formulation relate to the cell volume and cell face areas in the finite-volume 
formulation. For instance, because of equations (11), 


= Cell volume 


&r £y 

J' J’ J 


= x,y,z components of directed area of cell face perpendicular to 


£ = Constant coordinate plane 


which gives 





where are the unit vectors in the x,y,z directions, the subscript i refers to the cell- 

centered location (i, j, k), i- f- ^ corresponds to a cell interface location ( i -I- j, fc), and S is the 
directed area of the cell face normal to the £ coordinate. 


Numerical Algorithm 


Finite-volume solution procedure . Although references are made to differencing in the 
following discussion of the computational algorithm, it should be remembered that the actual 
operations are carried out by using the finite-volume method. In other words, the spatial 
derivatives shown in the equations are actually evaluated as a conservative flux balance across 
a cell. For example, 


OF 

d£ 


= «<(F) = F i+ , 


P H 


= (S-F) j+i — (S.F ) h 


(13) 


where 

F = Fl+Gj + Hfc 

By approximating the equations in this manner, one discretizes the equations for the 
numerical solution procedure. The interface fluxes S-F depicted in equation (13) are determined 
from the following characteristic-based method. 
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Interface flux . The convective and pressure terms at the cell interfaces are determined with 
the Roe type of upwind-biased flux-difference splitting of references 27 and 29. The method 
accounts for all different waves by which neighboring cells interact, including entropy and shear 
waves. In this method, the interface flux is written as the exact solution to an approximate 
Riemann problem 



1 ■ 
2 . 


F(Ql) + F(Q r ) - \A\(Q r - Q l ) 



(14) 


where A is A evaluated at Roe-averaged variables defined subsequently. Hence, 


A| = |A(Q)| 


and 


Q = Q(Ql,Qr) 


The variables Qi and Q# are the state variables to the left and right of the cell faces, and 


A = — = TAT -1 = T(A + + A - )T -1 

± = A±JA| 

2 

|A| = T|A|T -1 


(15) 


The diagonal matrix A is the matrix of eigenvalues of A, T is the matrix of right eigenvectors 
as columns, and T -1 the matrix of left eigenvectors as rows. They are all evaluated at Roe- 
averaged values such that 


F(Q* ) - F(Q L ) = A(Q* - Q X ) 


(16) 


is satisfied exactly. 


The third term in equation (14), |A|(Q# - Q L ), can be written as 
|A|(Qa - Ql) = |A|AQ 

0:4 

wa 4 + k x a 5 4- ag 
vct4 4- kya 5 4- a 7 
w)a 4 4- k z a 5 4- ag 

[ #«4 4- wag 4- wag + vaj + was - (^=t) a l J 


- < 


where 


a 1 = 


grad^ 


ia 


“ 2= 2p 


“ 3= 2^ 


grad e 


J 

grad e 


|u + c| (Ap + pcAu) 
\u - c\ (A p — pcAu) 


&4 = a l + a2 4- a3 
a 5 = c(a2 - <13) 

I gradC 


C*6 = 


a? = 


“8 = 


J 

grade 


J 

grade 


u\ (pAu — k x pAu ) 
S| (pAv — kypAu) 
u\ (pAw — k z pAu) 


The notation ~ denotes the following Roe-averaged evaluations 

P = \/PLPR 


~ _ u L + ur\J PrI PL 

1 + \ZprTpl 


~ _ n, + P /?/ pl 

1 + \TprFpl 


w L + WRsfpRfpL 

W = 


H = 


1 + \Jpr!pl 

Hl + H r^/prFpl 
1 + \FprTpl 


~2 / 1 \ TT U 2 + V 2 + W 2 

c 2 = h-l)H 


( 17 ) 


(18) 


( 19 ) 
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In equations (19), the magnitude of the cell interface directed area is 
cosines are 


{kxi ky,k z ) 


I grad £ | 


| gr | and the direction 

( 20 ) 


where the normalized contravariant velocity component normal to the cell interface is denoted 


u = k x u + k y v + k z w 


( 21 ) 


State-variable interpolations . The state-variable interpolations determine the resulting 
accuracy of the scheme. The state variables at the interface are constructed from the primitive 
variables q = [p, u, v, w,p] T . The interpolation 


(qi) i+ i = q* 
(q«) i+ i = q*+i ^ 


(22) 


corresponds to first-order fully upwind differencing. Higher order accuracy is given by the 
family of interpolations 


(qz,) i+ i = q» + \ [(1 - «)Vqj + (1 + «)Aqj] 
(q*)i+i = qi+l - \ [(! + «)Vqj+l + (1 - «)Sqj+l] 


(23) 


where Vq^ and Aq^ are, respectively, backward and forward differences of q. All the calculations 
i presented were made with k set to 1/3, a value which corresponds to a third-order accurate 

| spatial discretization. 

In order to ensure a monotone interpolation of q at the interfaces, the min mod flux limiter 
described in reference 11 

f Vq = min mod(Vq, 6Aq) 1 

_ \ ( 24 ) 
t Aq — min mod( Aq, bVq) J 


was used where 


\ min mod(x, y) = max {0, min[a: sign (y), by sign(x)]} sign(x) (25) 

and 

I 6=^ 

' 1 — K 


Time-differencing algorithm. The time-differencing algorithm used in the computational 
procedure is an approximately factored alternating-direction implicit scheme in delta form. 
Similar to the Beam and Warming technique of reference 30, equation (1) is factored and 
solved in the following steps: 
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1 s dF 

JAt + dQ 


i dG 

JAt + %Q 


AQ** = -R" 


AQ* = — AQ~ 


I 9H dH v 

JAt c l dQ dQ 


AQ - 7aJ aq * 


where AQ is the difference in Q between time levels n + 1 and n 

Q n+1 = Q n + AQ 

and 

R n = 6^F n + 6^ G n + <5 C (H - H„) n 


(26) 


(27) 

(28) 


The right-hand, or explicit, side of equations (26) is computed by evaluating the fluxes at 
the cell interfaces with the third-order accurate upwind-biased discretization. In the implicit, 
or left-hand side of the equation, the fluxes are evaluated with the first-order-accurate upwind 
discretization. 

The Jacobian matrices are evaluated explicitly. An example of their definition using the 
flux vector F is 

dF _ dF dQi ^ dF dQ R 

dQ dQ L dQ dQ R dQ [ ’ 

At a specific interface, dQi/dQ and dQ R /dQ determine the proper cell from which to obtain 
Q in computing each term in equation (29). For instance, for the i + % interface and using 
first-order-accurate upwind discretization, dQi/dQ is equal to 1 for a cell centered at i and is 
0 for all other locations. Futhermore, dQ R /dQ is 1 for a cell centered at i + 1 and 0 for all 
other locations. Thus, at the i + ^ interface, the first term in equation (29) becomes 


or 



9Ql 


4{^ + ?*-|akq r -q4 + ,} 


(30) 


dF. 


t+ 7 


dQ L 


1 

2 



m 

dQ L 


(Qr - Q l) + |A| i+ i 


(31) 


The second term in equation (31) is very small and is neglected. Approximating the third term 
gives 

dF ,1 , 

Similarly, 

9F i+ i ! 
dQ R ~ 2 ^ A * +1 ~ 

The Jacobian matrices for the flux at the i - j interface are determined following the same 
general procedure. 
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Substituting values similar to these spatially first-order-accurate upwind discretizations into 
equations (26) results in a block-tridiagonal system of equations to be solved in each of the 
three directions. In the i direction, this system of equations takes the form 

-A + (Qi-i, (S) f _ i ) 

+ tak + A+ ( Qi - (S M) - A_ ( Q - (S) h)] aq; * 

+ A" (Q i+ i,(S). +i ) AQfo = -R n (33) 

where 

A+ = i[A + |A|] 

A' = i (A - JA|] 

Each system of equations is solved with a full 5- by 5-block-tridiagonal inversion procedure. 

Boundary Conditions 

At the body, the velocity normal to the surface is set to zero and no-slip conditions are 
imposed to yield u = v = w = 0. An adiabatic wall condition is employed for temperature at 
the body. The pressure on the body surface is obtained by setting the normal pressure gradient 
equal to zero (a good assumption for viscous flow)and then extrapolating the pressure to the 
surface. 

The other boundaries of the computational domain are the inflow boundary, the far-field 
boundary, the outflow boundary at the downstream end of the domain, the axis ahead of the 
body, and the symmetry boundaries. Figure 1 illustrates these boundaries. At the inflow 
boundary, the treatment is based on the Riemann invariants for one-dimensional flow normal 
to the boundary. This method of treating a boundary is discussed thoroughly in references 31 
and 32. The same method of treatment is used at the far-field boundary where the flow at each 
individual point could be either into or out of the computational domain. At the downstream 
boundary, zeroth order extrapolation of all variables is used. Reflection conditions are used at 
the symmetry boundaries and on the axis. 
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Figure 1. Three-dimensional computational region and grid. 



Figure 2. Nonaxisymmetric computational body with solid jet exhaust plume simulator 
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Figure 3. Model mounted in wind tunnel. 
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(b) Sidewall geometry (centerline cross section). 






(c) Nozzle external cross-sectional geometry. 
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Figure 5. Concluded. 



Figure 6. Convergence properties and solution at Mach number of 0.80. Turbulent solution; R i = 19.5 x 10 6 ; 
coarse grid. 
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Figure 7. Continued. 
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Figure 7. Concluded. 








(d) Pressure distributions, side row. 




(f) Skin friction distributions, side row. 
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(d) Pressure distributions, side row. 




(f) Skin friction distributions, side row. 
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(e) Skin friction distributions, top row. 







Flow 



50 


Figure 12. Particle traces next to surface (computed oil flows). M 0 0 = 0.80; turbulent flow; R( = 19.5 x 10 6 ; 
fine grid. 




52 


Figure 13. Cross-flow velocity vectors. Looking upstream; M 0 0 = 0.80; turbulent flow; Rf = 19.5 x 10 6 ; fine 
grid. 
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Figure 13. Continued. 
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Figure 13. Continued. 
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Figure 13. Concluded. 
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